############################################################
#
#  RNPL program to solve 2 eqns resulting from the sourced 
#  3-D wave eqn.
# 
############################################################


#-----------------------------------------------------------
# Definition of memory size (only needed for Fortran)
#-----------------------------------------------------------
system parameter int memsiz := 250000000


#-----------------------------------------------------------
# Definition of parameters and associated default values. 
#-----------------------------------------------------------


#-----------------------------------------------------------
# The following parameters are used in the
# specification of the initial data. 
#-----------------------------------------------------------
parameter float A          := 0.1
parameter float deltax     := 0.1
parameter float deltay     := 0.1
parameter float deltaz     := 0.1
parameter float x0         := 0.0
parameter float y0         := 0.0
parameter float z0         := 0.0
parameter float epsKO      := 0.3 // Kreiss-Oliger dissipation parameter ( 0 <= epsKO < 1)
parameter float pi         := 3.14159265358979323846264338327950288

#-----------------------------------------------------------
# Note that "xmin" and "xmax" are special in that they are 
# also implicitly declared to be parameters via the 
# definition of the coordinate system below.
# Note to self: It turns out you can't declare these in terms
# of the quantities defined above. They must be numbers.
#-----------------------------------------------------------
parameter float xmin       := -1.0
parameter float xmax       :=  1.0
parameter float ymin       := -1.0
parameter float ymax       :=  1.0
parameter float zmin       := -1.0
parameter float zmax       :=  1.0


#-----------------------------------------------------------
# Definition of coordinate system, note that the first 
# coordinate is always assumed to be the time coordinate.
#-----------------------------------------------------------
rect coordinates t, x, y, z


#-----------------------------------------------------------
# Definition of finite-difference grids: [1:Nx] specifies
# the index range, {xmin:xmax} the coordinate range.
#-----------------------------------------------------------
uniform rect[x,y,z] grid g1 [1:Nx][1:Ny][1:Nz] {xmin:xmax}{ymin:ymax}{zmin:zmax}

#-----------------------------------------------------------
# Definition of grid functions: since a Crank-Nicholson 
# scheme is being used, the grid functions are defined at 
# "temporal offsets" 0 and 1, corresponding to "current" and 
# "advanced" time levels.  The directive {out_gf = 1} enables
# default output of the grid function (output can be disabled
# via {out_gf = 0}, by omitting the directive completely, or 
# via modification of the file .rnpl.attributes prior to 
# program invocation).
#-----------------------------------------------------------
float phi  on g1 at 0,1 {out_gf = 1}
float mu   on g1 at 0,1 {out_gf = 1}

# We define the metric (UU) and the source function:
float sqrtming on g1 at 0,1 {out_gf = 1}
float gtt on g1 at 0,1 {out_gf = 1}
float gtx on g1 at 0,1 {out_gf = 1}
float gty on g1 at 0,1 {out_gf = 1}
float gtz on g1 at 0,1 {out_gf = 1}
float gxx on g1 at 0,1 {out_gf = 1}
float gxy on g1 at 0,1 {out_gf = 1}
float gxz on g1 at 0,1 {out_gf = 1}
float gyy on g1 at 0,1 {out_gf = 1}
float gyz on g1 at 0,1 {out_gf = 1}
float gzz on g1 at 0,1 {out_gf = 1} 
float  FF on g1 at 0,1 {out_gf = 1}

#-----------------------------------------------------------
# FINITE DIFFERENCE OPERATOR DEFINITIONS
#-----------------------------------------------------------


#-----------------------------------------------------------
# Forward in time operator (for sanity check use)
#-----------------------------------------------------------
operator FORW(f,t) := <1>f[0][0][0]

#-----------------------------------------------------------
# Crank Nicholson time derivative operator (first forward 
# difference)
#-----------------------------------------------------------
operator DCN(f,t) := (<1>f[0][0][0] - <0>f[0][0][0]) / dt    

#-----------------------------------------------------------
# Forward time averaging operator
#-----------------------------------------------------------
operator AVE(f,t) := (<1>f[0][0][0] + <0>f[0][0][0]) / 2


# OPERATORS FOR X:
#-----------------------------------------------------------
# O(h^2) centred spatial derivative operator
#-----------------------------------------------------------
operator D0X(f,x) := (<0>f[1][0][0] - <0>f[-1][0][0]) / (2*dx)
operator DDCX(f,x) := (<0>f[1][0][0] - 2*<0>f[0][0][0] + <0>f[-1][0][0]) / (dx^2)


# OPERATORS FOR Y:
#-----------------------------------------------------------
# O(h^2) centred spatial derivative operator
#-----------------------------------------------------------
operator D0Y(f,y) := (<0>f[0][1][0] - <0>f[0][-1][0]) / (2*dy)
operator DDCY(f,y) := (<0>f[0][1][0] - 2*<0>f[0][0][0] + <0>f[0][-1][0]) / (dy^2)



# OPERATORS FOR Z:
#-----------------------------------------------------------
# O(h^2) centred spatial derivative operator
#-----------------------------------------------------------
operator D0Z(f,z) := (<0>f[0][0][1] - <0>f[0][0][-1]) / (2*dz)
operator DDCZ(f,z) := (<0>f[0][0][1] - 2*<0>f[0][0][0] + <0>f[0][0][-1]) / (dz^2)



#-----------------------------------------------------------
# Kreiss-Oliger dissipation operator
#-----------------------------------------------------------
operator DISS_KO(f,t) := - (epsKO/(16*dt))* (<0>f[-2][0][0] - 4*<0>f[-1][0][0] + 6*<0>f[0][0][0] - 4*<0>f[1][0][0] + <0>f[2][0][0] + 
                                             <0>f[0][-2][0] - 4*<0>f[0][-1][0] + 6*<0>f[0][0][0] - 4*<0>f[0][1][0] + <0>f[0][2][0] +
                                             <0>f[0][0][-2] - 4*<0>f[0][0][-1] + 6*<0>f[0][0][0] - 4*<0>f[0][0][1] + <0>f[0][0][2]
                                           )



#-----------------------------------------------------------
# DIFFERENCE EQUATION DEFINITIONS 
# Note the use of Dirichlet boundary conditions in that 
# phi and mu are defined to be 0 on the cube that encloses 
# the region.
#-----------------------------------------------------------

evaluate residual phi 
{ 
     [1:1][1:Ny][1:Nz]        := FORW(phi,t)  = 0;
     [1:Nx][1:1][1:Nz]        := FORW(phi,t)  = 0;
     [1:Nx][1:Ny][1:1]        := FORW(phi,t)  = 0;

     [2:2][2:Ny-1][2:Nz-1]    := DCN(phi,t)   = - AVE(  ( mu + gtx*D0X(phi,x) + gty*D0Y(phi,y) + gtz*D0Z(phi,z)   ) / gtt ,t);
     [2:Nx-1][2:2][2:Nz-1]    := DCN(phi,t)   = - AVE(  ( mu + gtx*D0X(phi,x) + gty*D0Y(phi,y) + gtz*D0Z(phi,z)   ) / gtt ,t);    
     [2:Nx-1][2:Ny-1][2:2]    := DCN(phi,t)   = - AVE(  ( mu + gtx*D0X(phi,x) + gty*D0Y(phi,y) + gtz*D0Z(phi,z)   ) / gtt ,t);  

     [3:Nx-2][3:Ny-2][3:Nz-2] := DCN(phi,t)   = - AVE(  ( mu + gtx*D0X(phi,x) + gty*D0Y(phi,y) + gtz*D0Z(phi,z)   ) / gtt ,t)  + DISS_KO(phi,t);

     [Nx-1:Nx-1][2:Ny-1][2:Nz-1]    := DCN(phi,t)   = - AVE(  ( mu + gtx*D0X(phi,x) + gty*D0Y(phi,y) + gtz*D0Z(phi,z)   ) / gtt ,t);
     [2:Nx-1][Ny-1:Ny-1][2:Nz-1]    := DCN(phi,t)   = - AVE(  ( mu + gtx*D0X(phi,x) + gty*D0Y(phi,y) + gtz*D0Z(phi,z)   ) / gtt ,t);     
     [2:Nx-1][2:Ny-1][Nz-1:Nz-1]    := DCN(phi,t)   = - AVE(  ( mu + gtx*D0X(phi,x) + gty*D0Y(phi,y) + gtz*D0Z(phi,z)   ) / gtt ,t);  

     [Nx:Nx][1:Ny][1:Nz]      := FORW(phi,t)  = 0;
     [1:Nx][Ny:Ny][1:Nz]      := FORW(phi,t)  = 0;
     [1:Nx][1:Ny][Nz:Nz]      := FORW(phi,t)  = 0;
}

evaluate residual mu 
{ 
     [1:1][1:Ny][1:Nz]        := FORW(mu,t)  = 0;
     [1:Nx][1:1][1:Nz]        := FORW(mu,t)  = 0;
     [1:Nx][1:Ny][1:1]        := FORW(mu,t)  = 0;

     [2:2][2:Ny-1][2:Nz-1]    := DCN(mu,t)  = 
                                               //T1:
                                               - AVE(mu/sqrtming,t)*DCN(sqrtming,t) + AVE(
                                               //T2:
                                               -FF
                                               //T3:
                                               +D0X(sqrtming*gxx,x)*D0X(phi,x)/sqrtming       
                                               +D0Y(sqrtming*gxy,y)*D0X(phi,x)/sqrtming
                                               +D0Z(sqrtming*gxz,z)*D0X(phi,x)/sqrtming

                                               +D0X(sqrtming*gxy,x)*D0Y(phi,y)/sqrtming       
                                               +D0Y(sqrtming*gyy,y)*D0Y(phi,y)/sqrtming
                                               +D0Z(sqrtming*gyz,z)*D0Y(phi,y)/sqrtming
      
                                               +D0X(sqrtming*gxz,x)*D0Z(phi,z)/sqrtming       
                                               +D0Y(sqrtming*gyz,y)*D0Z(phi,z)/sqrtming
                                               +D0Z(sqrtming*gzz,z)*D0Z(phi,z)/sqrtming
                                               //T4:
                                               +gxx*DDCX(phi,x) 
                                               +2*gxy*D0X(D0Y(phi,y),x)
                                               +2*gxz*D0X(D0Z(phi,z),x)

                                               +gyy*DDCY(phi,y)
                                               +2*gyz*D0Y(D0Z(phi,z),y)

                                               +gzz*DDCZ(phi,z)
                                               //T5:
                                               -D0X(sqrtming*gtx*mu/gtt,x)/sqrtming
                                               -D0Y(sqrtming*gty*mu/gtt,y)/sqrtming
                                               -D0Z(sqrtming*gtz*mu/gtt,z)/sqrtming
                                               //T6:
                                               -D0X(sqrtming*gtx*gtx/gtt,x)*D0X(phi,x)/sqrtming
                                               -D0Y(sqrtming*gtx*gty/gtt,y)*D0X(phi,x)/sqrtming
                                               -D0Z(sqrtming*gtx*gtz/gtt,z)*D0X(phi,x)/sqrtming

                                               -D0X(sqrtming*gty*gtx/gtt,x)*D0Y(phi,y)/sqrtming
                                               -D0Y(sqrtming*gty*gty/gtt,y)*D0Y(phi,y)/sqrtming
                                               -D0Z(sqrtming*gty*gtz/gtt,z)*D0Y(phi,y)/sqrtming

                                               -D0X(sqrtming*gtz*gtx/gtt,x)*D0Z(phi,z)/sqrtming
                                               -D0Y(sqrtming*gtz*gty/gtt,y)*D0Z(phi,z)/sqrtming
                                               -D0Z(sqrtming*gtz*gtz/gtt,z)*D0Z(phi,z)/sqrtming
                                               //T7:
                                               -gtx*gtx*DDCX(phi,x)/gtt 
                                               -2*gtx*gty*D0X(D0Y(phi,y),x)/gtt
                                               -2*gtx*gtz*D0X(D0Z(phi,z),x)/gtt

                                               -gty*gty*DDCY(phi,y)/gtt
                                               -2*gty*gtz*D0Y(D0Z(phi,z),y)/gtt  
 
                                               -gtz*gtz*D0Z(phi,z)/gtt
                                                                                       ,t);
     [2:Nx-1][2:2][2:Nz-1]    :=  DCN(mu,t) = 
                                               //T1:
                                               - AVE(mu/sqrtming,t)*DCN(sqrtming,t) + AVE(
                                               //T2:
                                               -FF
                                               //T3:
                                               +D0X(sqrtming*gxx,x)*D0X(phi,x)/sqrtming       
                                               +D0Y(sqrtming*gxy,y)*D0X(phi,x)/sqrtming
                                               +D0Z(sqrtming*gxz,z)*D0X(phi,x)/sqrtming

                                               +D0X(sqrtming*gxy,x)*D0Y(phi,y)/sqrtming       
                                               +D0Y(sqrtming*gyy,y)*D0Y(phi,y)/sqrtming
                                               +D0Z(sqrtming*gyz,z)*D0Y(phi,y)/sqrtming
      
                                               +D0X(sqrtming*gxz,x)*D0Z(phi,z)/sqrtming       
                                               +D0Y(sqrtming*gyz,y)*D0Z(phi,z)/sqrtming
                                               +D0Z(sqrtming*gzz,z)*D0Z(phi,z)/sqrtming
                                               //T4:
                                               +gxx*DDCX(phi,x) 
                                               +2*gxy*D0X(D0Y(phi,y),x)
                                               +2*gxz*D0X(D0Z(phi,z),x)

                                               +gyy*DDCY(phi,y)
                                               +2*gyz*D0Y(D0Z(phi,z),y)

                                               +gzz*DDCZ(phi,z)
                                               //T5:
                                               -D0X(sqrtming*gtx*mu/gtt,x)/sqrtming
                                               -D0Y(sqrtming*gty*mu/gtt,y)/sqrtming
                                               -D0Z(sqrtming*gtz*mu/gtt,z)/sqrtming
                                               //T6:
                                               -D0X(sqrtming*gtx*gtx/gtt,x)*D0X(phi,x)/sqrtming
                                               -D0Y(sqrtming*gtx*gty/gtt,y)*D0X(phi,x)/sqrtming
                                               -D0Z(sqrtming*gtx*gtz/gtt,z)*D0X(phi,x)/sqrtming

                                               -D0X(sqrtming*gty*gtx/gtt,x)*D0Y(phi,y)/sqrtming
                                               -D0Y(sqrtming*gty*gty/gtt,y)*D0Y(phi,y)/sqrtming
                                               -D0Z(sqrtming*gty*gtz/gtt,z)*D0Y(phi,y)/sqrtming

                                               -D0X(sqrtming*gtz*gtx/gtt,x)*D0Z(phi,z)/sqrtming
                                               -D0Y(sqrtming*gtz*gty/gtt,y)*D0Z(phi,z)/sqrtming
                                               -D0Z(sqrtming*gtz*gtz/gtt,z)*D0Z(phi,z)/sqrtming
                                               //T7:
                                               -gtx*gtx*DDCX(phi,x)/gtt 
                                               -2*gtx*gty*D0X(D0Y(phi,y),x)/gtt
                                               -2*gtx*gtz*D0X(D0Z(phi,z),x)/gtt

                                               -gty*gty*DDCY(phi,y)/gtt
                                               -2*gty*gtz*D0Y(D0Z(phi,z),y)/gtt  
 
                                               -gtz*gtz*D0Z(phi,z)/gtt
                                                                                       ,t);   
     [2:Nx-1][2:Ny-1][2:2]    := DCN(mu,t)  = 
                                               //T1:
                                               - AVE(mu/sqrtming,t)*DCN(sqrtming,t) + AVE(
                                               //T2:
                                               -FF
                                               //T3:
                                               +D0X(sqrtming*gxx,x)*D0X(phi,x)/sqrtming       
                                               +D0Y(sqrtming*gxy,y)*D0X(phi,x)/sqrtming
                                               +D0Z(sqrtming*gxz,z)*D0X(phi,x)/sqrtming

                                               +D0X(sqrtming*gxy,x)*D0Y(phi,y)/sqrtming       
                                               +D0Y(sqrtming*gyy,y)*D0Y(phi,y)/sqrtming
                                               +D0Z(sqrtming*gyz,z)*D0Y(phi,y)/sqrtming
      
                                               +D0X(sqrtming*gxz,x)*D0Z(phi,z)/sqrtming       
                                               +D0Y(sqrtming*gyz,y)*D0Z(phi,z)/sqrtming
                                               +D0Z(sqrtming*gzz,z)*D0Z(phi,z)/sqrtming
                                               //T4:
                                               +gxx*DDCX(phi,x) 
                                               +2*gxy*D0X(D0Y(phi,y),x)
                                               +2*gxz*D0X(D0Z(phi,z),x)

                                               +gyy*DDCY(phi,y)
                                               +2*gyz*D0Y(D0Z(phi,z),y)

                                               +gzz*DDCZ(phi,z)
                                               //T5:
                                               -D0X(sqrtming*gtx*mu/gtt,x)/sqrtming
                                               -D0Y(sqrtming*gty*mu/gtt,y)/sqrtming
                                               -D0Z(sqrtming*gtz*mu/gtt,z)/sqrtming
                                               //T6:
                                               -D0X(sqrtming*gtx*gtx/gtt,x)*D0X(phi,x)/sqrtming
                                               -D0Y(sqrtming*gtx*gty/gtt,y)*D0X(phi,x)/sqrtming
                                               -D0Z(sqrtming*gtx*gtz/gtt,z)*D0X(phi,x)/sqrtming

                                               -D0X(sqrtming*gty*gtx/gtt,x)*D0Y(phi,y)/sqrtming
                                               -D0Y(sqrtming*gty*gty/gtt,y)*D0Y(phi,y)/sqrtming
                                               -D0Z(sqrtming*gty*gtz/gtt,z)*D0Y(phi,y)/sqrtming

                                               -D0X(sqrtming*gtz*gtx/gtt,x)*D0Z(phi,z)/sqrtming
                                               -D0Y(sqrtming*gtz*gty/gtt,y)*D0Z(phi,z)/sqrtming
                                               -D0Z(sqrtming*gtz*gtz/gtt,z)*D0Z(phi,z)/sqrtming
                                               //T7:
                                               -gtx*gtx*DDCX(phi,x)/gtt 
                                               -2*gtx*gty*D0X(D0Y(phi,y),x)/gtt
                                               -2*gtx*gtz*D0X(D0Z(phi,z),x)/gtt

                                               -gty*gty*DDCY(phi,y)/gtt
                                               -2*gty*gtz*D0Y(D0Z(phi,z),y)/gtt  
 
                                               -gtz*gtz*D0Z(phi,z)/gtt
                                                                                       ,t);

     [3:Nx-2][3:Ny-2][3:Nz-2] :=  DCN(mu,t) - DISS_KO(mu,t) = 
                                               //T1:
                                               - AVE(mu/sqrtming,t)*DCN(sqrtming,t) + AVE(
                                               //T2:
                                               -FF
                                               //T3:
                                               +D0X(sqrtming*gxx,x)*D0X(phi,x)/sqrtming       
                                               +D0Y(sqrtming*gxy,y)*D0X(phi,x)/sqrtming
                                               +D0Z(sqrtming*gxz,z)*D0X(phi,x)/sqrtming

                                               +D0X(sqrtming*gxy,x)*D0Y(phi,y)/sqrtming       
                                               +D0Y(sqrtming*gyy,y)*D0Y(phi,y)/sqrtming
                                               +D0Z(sqrtming*gyz,z)*D0Y(phi,y)/sqrtming
      
                                               +D0X(sqrtming*gxz,x)*D0Z(phi,z)/sqrtming       
                                               +D0Y(sqrtming*gyz,y)*D0Z(phi,z)/sqrtming
                                               +D0Z(sqrtming*gzz,z)*D0Z(phi,z)/sqrtming
                                               //T4:
                                               +gxx*DDCX(phi,x) 
                                               +2*gxy*D0X(D0Y(phi,y),x)
                                               +2*gxz*D0X(D0Z(phi,z),x)

                                               +gyy*DDCY(phi,y)
                                               +2*gyz*D0Y(D0Z(phi,z),y)

                                               +gzz*DDCZ(phi,z)
                                               //T5:
                                               -D0X(sqrtming*gtx*mu/gtt,x)/sqrtming
                                               -D0Y(sqrtming*gty*mu/gtt,y)/sqrtming
                                               -D0Z(sqrtming*gtz*mu/gtt,z)/sqrtming
                                               //T6:
                                               -D0X(sqrtming*gtx*gtx/gtt,x)*D0X(phi,x)/sqrtming
                                               -D0Y(sqrtming*gtx*gty/gtt,y)*D0X(phi,x)/sqrtming
                                               -D0Z(sqrtming*gtx*gtz/gtt,z)*D0X(phi,x)/sqrtming

                                               -D0X(sqrtming*gty*gtx/gtt,x)*D0Y(phi,y)/sqrtming
                                               -D0Y(sqrtming*gty*gty/gtt,y)*D0Y(phi,y)/sqrtming
                                               -D0Z(sqrtming*gty*gtz/gtt,z)*D0Y(phi,y)/sqrtming

                                               -D0X(sqrtming*gtz*gtx/gtt,x)*D0Z(phi,z)/sqrtming
                                               -D0Y(sqrtming*gtz*gty/gtt,y)*D0Z(phi,z)/sqrtming
                                               -D0Z(sqrtming*gtz*gtz/gtt,z)*D0Z(phi,z)/sqrtming
                                               //T7:
                                               -gtx*gtx*DDCX(phi,x)/gtt 
                                               -2*gtx*gty*D0X(D0Y(phi,y),x)/gtt
                                               -2*gtx*gtz*D0X(D0Z(phi,z),x)/gtt

                                               -gty*gty*DDCY(phi,y)/gtt
                                               -2*gty*gtz*D0Y(D0Z(phi,z),y)/gtt  
 
                                               -gtz*gtz*D0Z(phi,z)/gtt
                                                                                       ,t);

     [Nx-1:Nx-1][2:Ny-1][2:Nz-1]    := DCN(mu,t)  = 
                                               //T1:
                                               - AVE(mu/sqrtming,t)*DCN(sqrtming,t) + AVE(
                                               //T2:
                                               -FF
                                               //T3:
                                               +D0X(sqrtming*gxx,x)*D0X(phi,x)/sqrtming       
                                               +D0Y(sqrtming*gxy,y)*D0X(phi,x)/sqrtming
                                               +D0Z(sqrtming*gxz,z)*D0X(phi,x)/sqrtming

                                               +D0X(sqrtming*gxy,x)*D0Y(phi,y)/sqrtming       
                                               +D0Y(sqrtming*gyy,y)*D0Y(phi,y)/sqrtming
                                               +D0Z(sqrtming*gyz,z)*D0Y(phi,y)/sqrtming
      
                                               +D0X(sqrtming*gxz,x)*D0Z(phi,z)/sqrtming       
                                               +D0Y(sqrtming*gyz,y)*D0Z(phi,z)/sqrtming
                                               +D0Z(sqrtming*gzz,z)*D0Z(phi,z)/sqrtming
                                               //T4:
                                               +gxx*DDCX(phi,x) 
                                               +2*gxy*D0X(D0Y(phi,y),x)
                                               +2*gxz*D0X(D0Z(phi,z),x)

                                               +gyy*DDCY(phi,y)
                                               +2*gyz*D0Y(D0Z(phi,z),y)

                                               +gzz*DDCZ(phi,z)
                                               //T5:
                                               -D0X(sqrtming*gtx*mu/gtt,x)/sqrtming
                                               -D0Y(sqrtming*gty*mu/gtt,y)/sqrtming
                                               -D0Z(sqrtming*gtz*mu/gtt,z)/sqrtming
                                               //T6:
                                               -D0X(sqrtming*gtx*gtx/gtt,x)*D0X(phi,x)/sqrtming
                                               -D0Y(sqrtming*gtx*gty/gtt,y)*D0X(phi,x)/sqrtming
                                               -D0Z(sqrtming*gtx*gtz/gtt,z)*D0X(phi,x)/sqrtming

                                               -D0X(sqrtming*gty*gtx/gtt,x)*D0Y(phi,y)/sqrtming
                                               -D0Y(sqrtming*gty*gty/gtt,y)*D0Y(phi,y)/sqrtming
                                               -D0Z(sqrtming*gty*gtz/gtt,z)*D0Y(phi,y)/sqrtming

                                               -D0X(sqrtming*gtz*gtx/gtt,x)*D0Z(phi,z)/sqrtming
                                               -D0Y(sqrtming*gtz*gty/gtt,y)*D0Z(phi,z)/sqrtming
                                               -D0Z(sqrtming*gtz*gtz/gtt,z)*D0Z(phi,z)/sqrtming
                                               //T7:
                                               -gtx*gtx*DDCX(phi,x)/gtt 
                                               -2*gtx*gty*D0X(D0Y(phi,y),x)/gtt
                                               -2*gtx*gtz*D0X(D0Z(phi,z),x)/gtt

                                               -gty*gty*DDCY(phi,y)/gtt
                                               -2*gty*gtz*D0Y(D0Z(phi,z),y)/gtt  
 
                                               -gtz*gtz*D0Z(phi,z)/gtt
                                                                                       ,t);
     [2:Nx-1][Ny-1:Ny-1][2:Nz-1]    :=  DCN(mu,t)  = 
                                               //T1:
                                               - AVE(mu/sqrtming,t)*DCN(sqrtming,t) + AVE(
                                               //T2:
                                               -FF
                                               //T3:
                                               +D0X(sqrtming*gxx,x)*D0X(phi,x)/sqrtming       
                                               +D0Y(sqrtming*gxy,y)*D0X(phi,x)/sqrtming
                                               +D0Z(sqrtming*gxz,z)*D0X(phi,x)/sqrtming

                                               +D0X(sqrtming*gxy,x)*D0Y(phi,y)/sqrtming       
                                               +D0Y(sqrtming*gyy,y)*D0Y(phi,y)/sqrtming
                                               +D0Z(sqrtming*gyz,z)*D0Y(phi,y)/sqrtming
      
                                               +D0X(sqrtming*gxz,x)*D0Z(phi,z)/sqrtming       
                                               +D0Y(sqrtming*gyz,y)*D0Z(phi,z)/sqrtming
                                               +D0Z(sqrtming*gzz,z)*D0Z(phi,z)/sqrtming
                                               //T4:
                                               +gxx*DDCX(phi,x) 
                                               +2*gxy*D0X(D0Y(phi,y),x)
                                               +2*gxz*D0X(D0Z(phi,z),x)

                                               +gyy*DDCY(phi,y)
                                               +2*gyz*D0Y(D0Z(phi,z),y)

                                               +gzz*DDCZ(phi,z)
                                               //T5:
                                               -D0X(sqrtming*gtx*mu/gtt,x)/sqrtming
                                               -D0Y(sqrtming*gty*mu/gtt,y)/sqrtming
                                               -D0Z(sqrtming*gtz*mu/gtt,z)/sqrtming
                                               //T6:
                                               -D0X(sqrtming*gtx*gtx/gtt,x)*D0X(phi,x)/sqrtming
                                               -D0Y(sqrtming*gtx*gty/gtt,y)*D0X(phi,x)/sqrtming
                                               -D0Z(sqrtming*gtx*gtz/gtt,z)*D0X(phi,x)/sqrtming

                                               -D0X(sqrtming*gty*gtx/gtt,x)*D0Y(phi,y)/sqrtming
                                               -D0Y(sqrtming*gty*gty/gtt,y)*D0Y(phi,y)/sqrtming
                                               -D0Z(sqrtming*gty*gtz/gtt,z)*D0Y(phi,y)/sqrtming

                                               -D0X(sqrtming*gtz*gtx/gtt,x)*D0Z(phi,z)/sqrtming
                                               -D0Y(sqrtming*gtz*gty/gtt,y)*D0Z(phi,z)/sqrtming
                                               -D0Z(sqrtming*gtz*gtz/gtt,z)*D0Z(phi,z)/sqrtming
                                               //T7:
                                               -gtx*gtx*DDCX(phi,x)/gtt 
                                               -2*gtx*gty*D0X(D0Y(phi,y),x)/gtt
                                               -2*gtx*gtz*D0X(D0Z(phi,z),x)/gtt

                                               -gty*gty*DDCY(phi,y)/gtt
                                               -2*gty*gtz*D0Y(D0Z(phi,z),y)/gtt  
 
                                               -gtz*gtz*D0Z(phi,z)/gtt
                                                                                       ,t);    
     [2:Nx-1][2:Ny-1][Nz-1:Nz-1]    := DCN(mu,t)  = 
                                               //T1:
                                               - AVE(mu/sqrtming,t)*DCN(sqrtming,t) + AVE(
                                               //T2:
                                               -FF
                                               //T3:
                                               +D0X(sqrtming*gxx,x)*D0X(phi,x)/sqrtming       
                                               +D0Y(sqrtming*gxy,y)*D0X(phi,x)/sqrtming
                                               +D0Z(sqrtming*gxz,z)*D0X(phi,x)/sqrtming

                                               +D0X(sqrtming*gxy,x)*D0Y(phi,y)/sqrtming       
                                               +D0Y(sqrtming*gyy,y)*D0Y(phi,y)/sqrtming
                                               +D0Z(sqrtming*gyz,z)*D0Y(phi,y)/sqrtming
      
                                               +D0X(sqrtming*gxz,x)*D0Z(phi,z)/sqrtming       
                                               +D0Y(sqrtming*gyz,y)*D0Z(phi,z)/sqrtming
                                               +D0Z(sqrtming*gzz,z)*D0Z(phi,z)/sqrtming
                                               //T4:
                                               +gxx*DDCX(phi,x) 
                                               +2*gxy*D0X(D0Y(phi,y),x)
                                               +2*gxz*D0X(D0Z(phi,z),x)

                                               +gyy*DDCY(phi,y)
                                               +2*gyz*D0Y(D0Z(phi,z),y)

                                               +gzz*DDCZ(phi,z)
                                               //T5:
                                               -D0X(sqrtming*gtx*mu/gtt,x)/sqrtming
                                               -D0Y(sqrtming*gty*mu/gtt,y)/sqrtming
                                               -D0Z(sqrtming*gtz*mu/gtt,z)/sqrtming
                                               //T6:
                                               -D0X(sqrtming*gtx*gtx/gtt,x)*D0X(phi,x)/sqrtming
                                               -D0Y(sqrtming*gtx*gty/gtt,y)*D0X(phi,x)/sqrtming
                                               -D0Z(sqrtming*gtx*gtz/gtt,z)*D0X(phi,x)/sqrtming

                                               -D0X(sqrtming*gty*gtx/gtt,x)*D0Y(phi,y)/sqrtming
                                               -D0Y(sqrtming*gty*gty/gtt,y)*D0Y(phi,y)/sqrtming
                                               -D0Z(sqrtming*gty*gtz/gtt,z)*D0Y(phi,y)/sqrtming

                                               -D0X(sqrtming*gtz*gtx/gtt,x)*D0Z(phi,z)/sqrtming
                                               -D0Y(sqrtming*gtz*gty/gtt,y)*D0Z(phi,z)/sqrtming
                                               -D0Z(sqrtming*gtz*gtz/gtt,z)*D0Z(phi,z)/sqrtming
                                               //T7:
                                               -gtx*gtx*DDCX(phi,x)/gtt 
                                               -2*gtx*gty*D0X(D0Y(phi,y),x)/gtt
                                               -2*gtx*gtz*D0X(D0Z(phi,z),x)/gtt

                                               -gty*gty*DDCY(phi,y)/gtt
                                               -2*gty*gtz*D0Y(D0Z(phi,z),y)/gtt  
 
                                               -gtz*gtz*D0Z(phi,z)/gtt
                                                                                       ,t);

     [Nx:Nx][1:Ny][1:Nz]      := FORW(mu,t)  = 0;
     [1:Nx][Ny:Ny][1:Nz]      := FORW(mu,t)  = 0;
     [1:Nx][1:Ny][Nz:Nz]      := FORW(mu,t)  = 0;
}


#-----------------------------------------------------------
# INITIALIZATION STATEMENTS 
#-----------------------------------------------------------

initialize phi 
{
     [1:Nx][1:Ny][1:Nz] := A* exp( -(x-x0)^2 / deltax^2 -(y-y0)^2 / deltay^2 -(z-z0)^2 / deltaz^2 )
}

initialize mu 
{
    [1:Nx][1:Ny][1:Nz] := - A* exp( -(x-x0)^2 / deltax^2 -(y-y0)^2 / deltay^2 -(z-z0)^2 / deltaz^2 )
}

initialize sqrtming { [1:Nx][1:Ny][1:Nz] := (8/pi^3)* ( cos(pi*x/2) )^2 * ( cos(pi*y/2) )^2 * ( cos(pi*z/2) )^2 }

initialize gtt { [1:Nx][1:Ny][1:Nz] :=-1 }
initialize gtx { [1:Nx][1:Ny][1:Nz] := 0 }
initialize gty { [1:Nx][1:Ny][1:Nz] := 0 }
initialize gtz { [1:Nx][1:Ny][1:Nz] := 0 }

initialize gxx { [1:Nx][1:Ny][1:Nz] := (4/pi^2) * ( cos(pi*x/2) )^4 }
initialize gxy { [1:Nx][1:Ny][1:Nz] := 0 }
initialize gxz { [1:Nx][1:Ny][1:Nz] := 0 }

initialize gyy { [1:Nx][1:Ny][1:Nz] := (4/pi^2) * ( cos(pi*y/2) )^4 }
initialize gyz { [1:Nx][1:Ny][1:Nz] := 0 }

initialize gzz { [1:Nx][1:Ny][1:Nz] := (4/pi^2) * ( cos(pi*z/2) )^4 }

initialize FF  { [1:Nx][1:Ny][1:Nz] := 0 }

auto initialize sqrtming, gtt, gtx, gty, gtz, gxx, gxy, gxz, gyy, gyz, gzz, FF, phi, mu


#-----------------------------------------------------------
# Definition of type of time stepping algorithm.  The use
# of "iterative" here, combined with the "auto update ..." 
# statement below, results in a scheme whereby the 
# residuals defined above are iteratively relaxed using 
# a point-wise Newton-Gauss-Seidel technique, until the 
# residual norms are below a certain threshold.
#-----------------------------------------------------------
looper iterative

#-----------------------------------------------------------
# The following statement directs RNPL to automatically 
# generate code to update the grid functions using the 
# residual definitions.
#-----------------------------------------------------------

metricandsource.inc update_metricandsource updates sqrtming, gtt, gtx, gty, gtz, gxx, gxy, gxz, gyy, gyz, gzz, FF
   header sqrtming, gtt, gtx, gty, gtz, gxx, gxy, gxz, gyy, gyz, gzz, FF, phi, mu, t, x, y, z

auto update phi, mu
